The development of new drug targets by the pharmaceutical industry is time-consuming and expensive because a large number of possible targets need to be tested before the molecule or compound with the desired properties is found or formulated. The same problem of a large number of possible variations affects the activities in the field of synthetic biology. In synthetic biology, biological entities are designed to perform a particular function, such as, for example, the development of biological nanomachines that may be used as programmed drug delivery systems. (See J. Panyam, V. Labhasetwar, Biodegradable nanoparticles for drug and gene delivery to cells and tissue, Advanced Drug Delivery Reviews, 55 (2003) 329-347.) As in drug discovery efforts, the formulation of a compound with desired properties is difficult due to the large variety of possible targets and the even larger context or system in which they must perform their function. Currently much of the work done to investigate the properties of these compounds is done in a wet-lab requiring many tedious and error prone experiments.
Development of chemical substances and nanomachinery, in addition to being time-consuming, can generate potentially dangerous intermediate substances. For example, a molecule used as a transport for a drug in a drug delivery system may, by its mere presence in the organism, stimulate the overproduction of some other protein. The overexpressed protein could act as a lethal toxin for the organism. Another possible complication is that the nanomachinery itself may mutate over time and either lose its original function or, worse, adversely interfere with the viability of the organism.
Another challenge facing the drug development activity is that, due to the cumbersome nature of experimental data collection, it is typical to limit experiments by narrowing the range of tested inputs and isolating the subsystem of interest. This limitation allows for the possibility that new drugs have unforeseen side-effects.
Moreover, current methods of obtaining data for biological processes are even more time-consuming than those associated with chemical processes, because the latter generally require laboratory experiments that lead to animal experiments and clinical trials. From these trials and experiments, data are obtained which, again, usually focus on a very narrow part of the biological system. Only after numerous costly trial-and-error clinical trials and constant redesigning of the clinical use of the drug to account for lessons learned from the most recent clinical trial, is a drug having adequate safety and efficacy finally realized. This process of clinical trial design and redesign, multiple clinical trials and, in some situations, multiple drug redesigns requires great expense of time and money. Even then, the effort may not produce a marketable drug. While conclusions may be drawn by assimilating experimental data and published information, it is difficult, if not impossible, to synthesize the relationships among all the available data and knowledge.
The various challenges faced by the aforementioned activities in chemical and biochemical research make it desirable to have software and methods for modeling, simulating, and analyzing biological processes in-silico rather than in-vitro or in-vivo. The goal of this approach is to provide a more comprehensive view of these biological systems prior to costly experiments and to clinical trials thereby reducing the search space for drug targets and useful nanoparticles.
Dynamic systems, such as biological processes and chemical reactions can be modeled as sets of differential, difference, algebraic, and/or recursive equations. At any given instant of time, these equations may be viewed as relationships between the system's output response (“outputs”), the system's input stimuli (“inputs”) at that time, the current state of the system, the system parameters, and time. The state of the system may be thought of as a numerical representation of the dynamically changing configuration of the system. For instance, in a physical system modeling a simple pendulum, the state may be viewed as the current position and velocity of the pendulum. Similarly, a signal-processing system that filters a signal would maintain a set of previous inputs as the state. The system parameters are the numerical representation of the static (unchanging) configuration of the system and may be viewed as constant coefficients in the system's equations. For the pendulum example, a parameter is the length of pendulum and for the filter example; a parameter is the values of the filter taps.
Inherent in four of the classes of systems (ODE, difference equations, algebraic equations and composite) is the notion of system sample time. The sample-time is the time interval at which the inputs, state, or outputs (collectively referred to as the results) of the system are traced as time progresses. Based on sample times, a system can be described as a discrete-time system, continuous-time system and hybrid system. As noted above, stochastic systems occur at a random time determined by a reaction-specific operative probability distribution.
A discrete-time system is a system in which the evolution of the system results is tracked at finite intervals of time. In the limit as the interval approaches zero, the discrete-time system becomes a continuous-time system. The intervals of time may be periodic or non-periodic. Sometimes, non-periodic rate systems, such as stochastic systems, are referred to as non-uniform rate systems meaning that there is no periodic rate at which the response can be tracked. A continuous-time system is a system in which the evolutions of the system results are continuously changing. Continuous-time signals change during numerical integration. An example of a continuous-time system is one described by an ODE. There can also be algebraic or composite continuous-time systems. A hybrid system is a system with both discrete-time and continuous-time elements.
If a system has only one sample time, it is said to be single-rate. If a system has multiple sample times, it is said to be multi-rate. Multi-rate systems can be evaluated (executed) using either a single-tasking form of execution or a multi-tasking form of execution. When multi-tasking execution is used, it conforms to rate monotonic scheduling principals as defined by Liu, C. L., and LAYLAND, J. W. Scheduling Algorithms for Multiprogramming in a Hard-Real-Time Environment. ACM 20, 1 (January 1973), 46-61. Systems may also be categorized by the type of numerical integration solver being used. A fixed-step system is one that uses a fixed-step solver. Fixed-step solvers typically use explicit methods to compute the next continuous state at fixed periodic intervals of time. A variable-step system is one that is using a variable-step solver. A variable-step solver can use either implicit or explicit methods to compute the next continuous state at non-periodic intervals of time. Generally, variable-step solvers use a form of error control to adjust the interval size such that the desired error tolerances are achieved.
For biological process and chemical reaction models, stochastic solvers may be useful because stochastic reactions occur at a random time based on a reaction-specific probability distribution and hence they do not neatly fit either a fixed-step type of solver or a continuous-time solver. The stochastic solvers may use the exact stochastic simulation algorithm (SSA).
The exact SSA numerically simulates the time evolution of a given chemical system. In the SSA technique, reaction events given selected probabilities of occurring, and the events which occur change the probabilities of subsequent events. The algorithm determines, for a system in a given state, the next reaction to occur and the time that the next reaction occurs using probability. The algorithm is based on a quantity P(t, u), which is the probability that a reaction u will occur at the time interval t. The probabilities are based on the classical rate coefficients (k), the volume of the container, which can be a cell, a partition of a cell, a compartment of the cell, such as the nucleus or other organelles, or other container, and the concentration of reactants in a given reaction. Once a time and reaction have been computed, the method carries out the reaction, i.e., it updates the state of the system to reflect the transformation of reactants into products, then increments the time by t and determines another reaction to occur and when the reaction will occur. The SSA technique is described in detail in the article: Gillespie, D. T. 1977, Exact Stochastic Simulation of Coupled Chemical Reactions, Journal of Physical Chemistry, vol. 81, pp. 2340-2361.
Since the SSA simulates every reaction event, the simulation result of the SSA is accurate but it is too slow for practical simulation of the chemical or biological reaction systems. The tau-leaping algorithms have been proposed to accelerate the SSA by leaping over sequences of non-critical reactions that occur in a time interval, tau (τ). In the tau-leaping algorithm, the size of the time interval is taken to encompass more than one reaction. Since the tau-leaping algorithms accelerate the SSA by sacrificing the accuracy of the simulation result, it is important to select the size of the time interval properly in order to balance the acceleration and accuracy of the simulation.
A recent article has proposed a method for determining the size of the time interval in an explicit tau-leaping algorithm, which is used to simulate non-stiff systems. See the article: D. T. Gillespie and L. R Petzold, Improved Leap-Size Selection for Accelerated Stochastic Simulation, Journal of Physical Chemistry, vol. 119 (2003), pp. 8229. The explicit tau-leaping algorithm, however, does not produce good results when applied to a stiff system. Therefore, an implicit tau-leaping algorithm has been proposed to address the stiffness of the system. See the article: T. Rathinam, L. R. Petzold, and D. T. Gillespie, Stiffness in Stochastic Chemically Reacting System: The Implicit Tau-Leaping Method, Journal of Physical Chemistry, vol. 119 (2003), pp. 11784-94. The article, however, does not provide a method for selecting the size of the time interval properly in the implicit tau-leaping algorithm. Automatic selection of the size of the time interval based on the model and the user input is necessary in a general purpose simulation tool.